#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
    Test de Primalidad Probabilístico Miller-Rabin
    Copyright (C) 2012 Darío López Padial.
    @bukosabino

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>
"""

################################################################################
#
#	Ejecución:
#		$ ./primos.py
#
#	Objetivo: 
#	El programa decidirá si un número es probable primo o no.
#	Podremos indicar la probabilidad de testeo para la prueba o el número de
#	testeos que deseamos. La respuesta será "no es primo" o "es probable primo".
#
################################################################################

import sys
import random

# Teorema de Fermat para calcular a^m (mod n)
def fermat(a, m, n):
	b = 1
	while(m!=1):
		if m%2 == 1:
			b = long((a*b) % n)
		a = pow(a, 2) % n
		m = m / 2
	b = long((a*b) % n)
	return b

# Test de primalidad de Miller-Rabin
def miller_rabin():
    num = random.randint(2, p-2)
    a = long(fermat(num, s, p))
    if a != 1 and a != (p-1):
        j=1
        while j<u:
            a = long(fermat(a, 2, p))
            if a == (p-1):
                return True
            if a == 1:
                return False
            j+=1
        return False
    else:
        return True	

	
def salida(primo, probabilidad):
    if primo == True:
	    if probabilidad == 100:
	    	probabilidad = float( float( pow(4, n) - 1 ) / pow(4, n) )
	    print "El número introducido es probablemente primo, con una probabilidad de", (probabilidad * 100)
    else:
    	print "El número introducido no es primo."
    sys.exit()
    

print u'Test probabilístico de primalidad, de Miller-Rabin: '
booleano = True
while(booleano):
	gen = raw_input('¿Desea generar aleatoriamente el "número primo" a analizar?¿s/n?')
	if gen == 's':
		booleano = False
		num = long(input('Señala el número de cifras que desea para el número :'))
		p = long(random.randint(pow(10, num-1), pow(10, num)))
		if p % 2 == 0:
			p += 1
		print u"El número escogido es", p
	elif gen == 'n':
		booleano = False
		p = long(input('Introduce el valor del "número primo" a analizar: '))
		if p % 2 == 0:
			print u"El número introducido es par."
			sys.exit()

probabilidad = 100
booleano = True
while(booleano):
	opc = raw_input('¿Desea establecer una probabilidad mínima con la que probar la primalidad?¿s/n?')
	if opc == 's':
		probabilidad = input('Indique numéricamente del 0 al 100 la probabilidad de seguridad que desea: ')
		n = 50 #valor grande, para no pisarnos.
		booleano = False
	elif opc == 'n':
		booleano = False
		n = input('Introduce un número natural, este cuantificará el número de test que se harán para comprobar la primalidad: ')

s = long(p - 1)
u = 0L
while (s%2==0):
	u += 1
	s = s/2
for i in range(n):
    prob = float( float( pow(4, i) - 1 ) / pow(4, i) )
    if probabilidad < prob*100:
        probabilidad = prob
        break
    if miller_rabin() == False:
        salida(False, probabilidad)
salida(True, probabilidad)
